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The problem of plane-wave diffraction on semi-infinite orthorhombic electromagnetic (photonic) 
crystals of general kind is considered. Boundary conditions are obtained in the form of infinite 
system of equations relating amplitudes of incident wave, eigenmodes excited in the crystal and 
scattered spatial harmonics. Generalized Ewald-Oseen extinction principle is formulated on the 
base of deduced boundary conditions. The knowledge of properties of infinite crystal's eigenmodes 
provides option to solve the diffraction problem for the corresponding semi-infinite crystal numer- 
ically. In the case when the crystal is formed by small inclusions which can be treated as point 
dipolar scatterers with fixed direction the problem admits complete rigorous analytical solution. 
The amplitudes of excited modes and scattered spatial harmonics are expressed in terms of the 
wave vectors of the infinite crystal by closed-form analytical formulae. The result is applied for 
study of reflection properties of metamaterial formed by cubic lattice of split-ring resonators. 
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I. INTRODUCTION 

Electromagnetic crystals are artificial periodical struc- 
tures operating at the wavelengths comparable with their 
periods [1-3]. At the optical frequencies such structures 
are called as photonic crystals. The inherent feature of 
these materials is the existence of frequency bands where 
the crystal does not support propagating waves. The 
band gaps are caused by spatial resonances of the crystal 
lattice and strongly depend on the direction of propaga- 
tion. It means that electromagnetic crystals are media 
with spatial dispersion [4-6]. The material parameters: 
permittivity and permeability for such materials, if they 
can be introduced at all, depend on the wave vector as 
well as on the frequency. Notice, that the homogenization 
approach is not the most convenient way for the descrip- 
tion of electromagnetic crystals even at low frequencies. 
It often requires introduction of additional boundary con- 
ditions in order to describe boundary problems correctly, 
and this involves related complexities. The photonic and 
electromagnetic crystals are usually studied with the help 
of numerical methods [1-3] . Analytical models exist only 
for a very narrow class of the crystals. Some types of 
the crystals can be studied analytically under a certain 
approximation, but the strict analytical solution for a 
photonic crystal is an exception. 

The goal of the present paper is to demonstrate how 
boundary problems for electromagnetic crystals can be 
effectively studied using analytical methods. The paper 
is separated into the two parts. In the first part the 
boundary conditions for electromagnetic crystals of gen- 
eral kind are deduced in the form of infinite system of 



equations relating amplitudes of incident wave, excited 
eigenmodes of the crystal and scattered spatial harmon- 
ics. This system can be interpreted as generalization of 
well-known Ewald-Oseen extinction principle [7-9] which 
states that the polarization of dielectric is distributed 
so that it cancels out the incident wave and produces 
the propagating wave. For the electromagnetic crystals, 
inherently periodic structures, the generalized Ewald- 
Oseen principle states that the polarization of dielectric 
is distributed so that it cancels out the incident wave as 
well as all spatial harmonics associated with periodicity 
of the boundary. This principle expressed in the form of 
infinite system of boundary conditions provides opportu- 
nity to solve the boundary problem for semi-infinite crys- 
tal of certain kind numerically if the eigenmode problem 
for corresponding infinite crystal is already solved. In 
the second part of the paper the proposed approach is 
applied for the case of electromagnetic crystals formed 
by small inclusions which can be treated as point dipolar 
scatterers with fixed direction. In this case the system 
of boundary conditions admits complete rigorous ana- 
lytical solution. The amplitudes of excited eigenmodes 
and scattered spatial harmonics are expressed in terms 
of wavevectors of eigenmodes using closed-form analyt- 
ical formulae. These results are unique extension and 
generalization of known Mahan-Obermair theory [10] for 
the case then period of the crystal is compared with wave- 
length. At the end of the paper it is demonstrated how 
reflection from the semi-infinite cubic lattice of resonant 
scatterers (split-ring resonators) can be modeled in the 
regime of strong spatial dispersion observed in such crys- 
tals [11]. 



II. PROOF OF GENERALIZED EWALD-OSEEN 
EXTINCTION PRINCIPLE 

In this section we provide proof of generalized Ewald- 
Oseen extinction principle for arbitrary semi-infinite elec- 
tromagnetic crystal with orthorhombic elementary cell. 
First, let us consider an infinite orthorhombic electro- 
magnetic crystal with geometry schematically presented 
in Fig. 1 and characterized by three-periodical permit- 
tivity distribution: 



e(r) = e(r + an + hs + cl). 



(1) 



In this expression and in the further text the two lines 
over a quantity designate that the quantity is dyadic (ten- 
sor of second rank in three-dimensional space) . It means 
that we consider of the most general kind of electromag- 
netic crystals formed by dielectrics. 
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FIG. 1: Geometry of an infinite electromagnetic crystal 

In this paper we are using local field approach, an un- 
conventional method for description of fields inside di- 
electrics. We will operate with local parameters like po- 
larization density P and local electrical field Ei oc ., but 
not with average electric field E and displacement D as 
usual. The similar approach was used in [9] for rigorous 
derivation of Ewald-Oseen extinction theorem and in [6] . 
The dielectric can be treated as a very dense cubic lat- 
tice of point scatterers with ceratin local polarizability. 
In this formulation the dielectric permittivity e(r) has to 
be replaced (see Figure 2 for illustration) by the local 
polarizability a(r) relating the bulk polarization density 
P(r) with the local electric field Ei oc (r): 



P(r) = a(r)Eioc(r). 



The expression for local polarizability in terms of di- 
electric permittivity has the following form: 



a (r) 



e(r) - e Q I 



J/(3eo) 



(3) 



where £o is permittivity of free space, / is unit dyadic. 
This expression follows from the Lorentz-Lorenz formula 
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FIG. 2: Illustration for replacement of dielectric permittivity 
by local polarizability 



[9] 



(4) 



Ei oc (r) = E(r)+P(r)/(3e ) 
and material equation 

D(r)=e E(r)+P(r)=f(r)E(r). (5) 

A. Dispersion equation 

Following local field approach one can write down dis- 
persion equation for the crystal under consideration in 
the next integral form: 

P(r) = 5(r) / <5 3 (r - r', q)P(r')dr', VreV, (6) 



where V = V(a, b, c) is volume of the elementary lattice 
cell, G3(r, q) is the lattice dyadic Green's function: 



bs+q z cl) 



G 3 (r, q) = G ( r - an - bs - c /)e^ ( ^ Qn+ ^ 



. (7) 

which takes into account cell-to-cell polarization distri- 
bution determined by wave vector q = (q x , q y , q z ) T : 



( 2 ) 55 



P(r + an + bs + cl) = V{Y)e- j{ ^ an+q y hs+q ^ cl \ 
G(y) is dyadic Green's function of free space: 

<5(r) = (fc 2 T+ VV) 



47T£o?" 



(8) 



(9) 



n, s,l are integer indices, k is wave number of free space. 
The integral in (6) is singular if the point corresponding 
to vector r is located inside of some polarized dielectric. 
It has to be evaluated in the meaning of principal value 
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by excluding small spherical region around the singular 
point and tending the radius of this region to zero [12]. 

The dispersion equation (6) relates distribution of po- 
larization density P(r) and wave vector q correspond- 
ing to the eigenmodes of the electromagnetic crystal. If 
the distribution of average electric field E(r) of a crystal 
eigenmode is known then the polarization density P(r) 
can be found directly using material equation (5): 

P(r) = P(r)- £o ]E(r). (10) 

The reverse operation is possible only for space regions 
filled by dielectric with e(r) ^ eq. The distribution of 
electric field in free space regions, if required, have to be 
calculated using the next integral representation 

E(r) = J 5 3 (r - r', q)P(r')dr'. (11) 

v 

For our proof of generalized Ewald-Oseen extinction 
principle we have to transform dispersion equation (6) 
into the form corresponding to summation by layers in 
x-direction. The expression for lattice dyadic Green's 

function Gs(r, q) (7) can be rewritten using summation 
over planes in the form: 

S 3 (r,q)= J2 %(r~an)e-^ an , (12) 

71— — OO 

where G2(r) is the grid dyadic Green's function: 

<5 2 (r) =J2^( r - hs - cl ) e ~ KqvbS+9zCl) ■ ( 13 ) 

s,l 

Applying Poisson summation formula by both indices s 
and I one can express the grid dyadic Green's function in 
terms of the spatial Floquet harmonics. This expansion 
is also called as spectral representation: 

^(r)^^tf l(a; " a) e^ (k - <X, - r) , (W) 

s,l 

where 

= ^J-^tl X Kl X ^ Kl = (±*?,!>*a W .*f) T » 

k y s = % + — , kf = q z + — , k x s l = \J k 2 - {hi) 2 - (kf ) 2 . 

The square root in the expression for fc S) ; should be cho- 
sen so that Im(^A) < 0. The sign ± corresponds to half 
spaces x > a and x < a respectively. 

Using (12) the dispersion equation (6) can be rewritten 
in the following form which will be used later on: 

P(r) = E(r) / ^(r - r' - an)P(r') e - j ^ an dr' . 



B. Semi-infinite crystal 

Now let us consider a semi-infinite crystal (half space 
x > a, see Figure 3) excited by a plane electromagnetic 
wave with wave vector k = (k x , k y , k z ) T coming from free 
space: 

E inc (r) =E inc e-^ k - r >. (16) 

The origin of our coordinate system is intentionally 
shifted by one period into the free space since it simpli- 
fies rather cumbersome calculations which are presented 
below and causes exponential convergence of series in the 
final expressions. 




FIG. 3: Geometry of an semi-infinite electromagnetic crystal 

Due to the periodicity of the semi-infinite structure 
along y— and z— axes the distribution of the excited po- 
larization along these directions is determined by the 
phase of the incident wave: 

P(r + am + bs + cl) = P(r + am)e^ (fc " b ™ + ^ c/) , (17) 

for any r € V and m > 1. 

The electric field produced by the polarized semi- 
infinite crystal has the form 

+oo . 

E scat (r) =J2 ^ r - r ' - a ") p ( r ' + a ™) dr '- ( 18 ) 
»=iv 

The total local electric field is sum of incident and scat- 
tered (produced by polarization of crystal) fields. Follow- 
ing the local field approach we can write: 

P(r) = 5(r) [E inc (r) + E scat (r)] . (19) 

Combining (18) and (19) we obtain an integral equa- 
tion for the polarization in the semi-infinite crystal ex- 
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cited by an incident wave: 



P(r) = a(r) 



+ 00 



E inc (r) + E / G ^ r - r ' - a «) p ( r ' + a «) rfr ' 

(20) 



Further, substituting (14) into (24), changing the sum- 
mation order and evaluating the sum of the geometrical 
progression by index n we get: 



7;jP i (rV (k - rr ' , rfr' 



Now let us suppose that the dispersion equation (6) 
is solved under condition that wave vector has the form 
q = (q x ,k y ,k z ) T with unknown a;— component q x , and 
the set of eigenmodes {(Pj(r), <j4^)}, characterized by x- 

(i) 

component of wave vector q x and polarization distri- 
bution P»(r), is found. In addition, we include into this 
set only the eigenmodes which either transfer energy into 
half-space x > a [dq x ^ /du> > 0] or decay in x-direction 
im < 0]. In such a case the polarization of the semi- 
infinite crystal excited by the incident wave (16) can be 
expanded by the eigenmodes of infinite crystal as follows: 



EE* 



-j(k+,.(r+am)) 



Ei, 



-jk(r-fam) 



(25) 



C. Generalized Ewald-Oseen extinction principle 

The left part of equation (25) represents an expansion 
of the right part into a spatial spectrum of Floquet har- 
monics. The right part represents an incident spectrum 
of Floquet harmonics containing only the single incident 
plane wave (16) with k = kj . Equating coefficients in 



E. m pi<anu w<avu I iu I Willi r*. — i^-nn- -L^M liciliiii; 

A^(r)e-^-, Vr e V,m > 1. (21) ^ ^ partg q{ (2 ° 5) ° we " obtain ? 



Substituting (21) into (20) and using (16) we obtain: 



t«E*- 



/P i (r')e j(k M- r ' ) rfr / 



1 _ e J'(9* ) -feJ,i)o 



_ / E inc , (a,l) = (0,0) 
0, (s,0^(0,0) 



J2AiPi(r)e- jq * )am = a(r) 



E inc e^ k(r+am) (22) 



+00 „ 

+EE* / 



G 2 (r - r' - a(n - ro))Pi(r')e-«- a "dr' 



Splitting series in the dispersion equation (15) one can 
derive the following auxiliary relation: 

P i (r)e-^ i)aro = 



S(r) J2 / G 2 (r - r' - a(n - m))P 4 (r')e- J ^ >a "dr' 

(23) 



n— 1 1 



+a(r) ^ / G 2 (r-r' -a(n-m))P 4 (r')e^'^ a "dr'. 



Substituting (23) into (22) and following the fact that 
dct{a(r)} ^Owe obtain: 

E* E / G 2 (r-r' -a(n-m))P 4 (r')e^ >a "dr' 



(26) 

The values in the right side of (26) are the amplitudes 
of the incident spatial harmonics (all harmonics except 
fundamental one have zero amplitudes), and the series 
in the left side are the amplitudes of the spatial har- 
monics produced by the whole semi-infinite crystal po- 
larization in order to cancel these incident harmonics. It 
means that equation (26) represents the generalization of 
Ewald-Oseen extinction principle (see [7-9] for classical 
formulation in the case of dielectrics): the polarization in 
a semi-infinite electromagnetic crystal excited by a plane 
wave is distributed in such a way that it cancels the inci- 
dent wave together with all high-order spatial harmonics 
associated with periodicity of the boundary (even if they 
have zero amplitudes as in the present case). The addi- 
tional words related to high-order Floquet harmonics is 
the main and principal difference of Ewald-Oseen extinc- 
tion principle formulation for electromagnetic crystals as 
compared to the classical case of isotropic dielectrics. 

Substitution of (21) and (14) into (18) allows to express 
the scattered field in the half space x < a in terms of 
spatial Floquet harmonics: 



E s 



E 



E scat = Is, I E* 



fPi(rV (k M- r 'W 



(27) 



(28) 



P r p -jk(r+am) 

J_J inc O 



(24) 



Note, that the formula for the amplitudes of scattered 
Floquet harmonics (28) contains series which have the 
same form as (26) and differs only by the sign of the 
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x— components of wave vectors k s l — (±k* ; , k y , k z ) T 
corresponding to the spatial harmonics propagating into 
the half spaces x < a and x > a, respectively. 

If the eigenmodes of the crystal {qi\ Pj( r )} are known 
then one can solve the system of linear equations (26) 
and find amplitudes of excited eigenmodes {Ai}. With 
use of these amplitudes the scattered field can be found 
by (28). This provides a new numerical method which 
allows to solve problem of plane-wave diffraction by a 
semi-infinite electromagnetic crystal using knowledge of 
eigenmodes of the infinite crystal. This fact is very im- 
portant since at the moment the reflection and dispersion 
problems for electromagnetic crystals are usually solved 
by separate numerical approaches. The expressions (26) 
and (28) create a link between these two problems and 
show how results of dispersion studies can be used in or- 
der to describe reflection properties of electromagnetic 
crystals. 

Until this point we considered only one incident wave 
with wave vector k = kj , but from (25) it is clear that 
we could consider also other incident spatial harmonics 
with wave vectors kj~ ; and get the similar results as (26) , 
but the nonzero terms at the right side of equation would 
correspond to the respective incident spatial harmonic. 
Using principle of superposition we obtain that if the 
semi-infinite crystal is excited by whole spectrum of in- 
cident spatial harmonics with amplitudes E^ c and wave 
vectors kj~ ; then the following system of linear equations 
is valid: 



/P J (r')e j(k ^- r ' ) ( ir' 



1 - e 



j(?£ -fef ,)<« 



= Ef 



(29) 



to have negative polarization density with respect to di- 
electrics with e > Eq. This can be simply explained since 
free space with respect to these dielectrics is like real ma- 
terials with £ < £o with respect to free space: they indeed 
have negative polarization density. 

The meaning of polarization density becomes relative 
automatically when the replacement of e to e is made. 
In order to avoid the use of this ambiguous polarization 
in the final formulae it is possible to express the polariza- 
tion density in terms of the average field using (5). The 
resulting expressions provide complete set of boundary 
conditions for interface between semi-infinite electromag- 
netic crystal and isotropic dielectric: 

' V fee"*."..' + KLe- j <' r ) ,x<a 
E(r) = { ^ V ' 
1 J2A l E l e-^ r , x>a 



J[t-£ /]E J (r')e-'' ( <'- r ' ) dr' 

I _ gj(9i - fe s,i)° 

/|f-eoT|Ei(r , )e J '( k ."'- p 'W 



(30) 



K c , (31) 



I _ gjfe +K,i) a 
where we use following notations 



Ef 



(32) 



ls,l = 



2bcek 



8,1 



[k± X [k± x/]], k± ={±ki l ,ky,ktf 



The scattered field is given by (28) as in the case of one 
incident wave. The equation (29) represents cancelation 
of all incident spatial spectrum by induced polarization 
of the crystal in accordance to the formulated above gen- 
eralized Ewald-Oseen extinction principle. 



D. Formulation of boundary conditions 

In the previous section we have proved generalized 
Ewald-Oseen extinction principle for the case of semi- 
infinite electromagnetic crystal described by certain pe- 
riodic permittivity distribution e(r) excited by a plane 
wave coming from free space. In order to extend this 
theory for the case when incident wave comes from ho- 
mogeneous isotropic dielectric with permittivity e it is 
enough to change eo in a U formulae to e. Physically it 
means that we have to consider the polarization of the 
crystal with respect to the host material with permittiv- 
ity £, but not free space. In the model of dense cubic 
lattice of point dipoles it means that the lattice is lo- 
cated inside of this host material. This approach is very 
unusual since it can lead to results which are strange 
from first point of view. For example, free space happen 



y _ 2ws 
k s — c ly~\ i H 



2ttZ 



\J k 2 - (fcf ) 2 



(fcf) 2 , 



k is wave number in the dielectric with permittivity e, 
1i = (if \k y ,k z ) T and the square root in the expression 
for k Sy i should be chosen so that Im(\A) < 0. 

The expressions (31) and (32) relate amplitudes of in- 
cident E^' c and scattered E^ z at spatial harmonics corre- 
sponding to tangential wave vector k t = (k y ,k z ) T and 
periodicity of the boundary (rectangular lattice with pe- 
riods b and c) with amplitudes Ai of eigenmodes (Ej, qf ) 
excited in the semi-infinite crystal. 

The presented set of boundary conditions is complete: 
these equations are enough to determine amplitudes of 
excited eigenmodes and scattered spatial harmonics if the 
eigenmodes {(E i7 qf )} of infinite crystal corresponding to 
tangential wave vector k t are known. But this set is not 
unique. One can immediately suggest to use classical 
boundary conditions (continuous tangential component 
of electric field and normal component of electric dis- 
placement at any point of the boundary x = a) which 
being expanded into Fourier series will also grant a com- 
plete set of linear equations relating amplitudes of inci- 
dent, scattered and excited modes. 
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The advantage of equations (31) and (32) as compared 
to any other boundary conditions is such that they have 
very special form which can admit analytical solution. 
We will demonstrate it in the next section for the special 
case of electromagnetic crystals formed by small scat- 
terers which can be treated as point dipolc with fixed 
orientation. But this is not the only case when an analyt- 
ical solution of (31) and (32) can be obtained. Recently, 
M. Silveirinha [13] demonstrated that the method pro- 
posed by ourselves can be successfully applied for stud- 
ies of reflection from semi-infinite wire medium, material 
with strong low-frequency spatial dispersion [14]. Unfor- 
tunately, we can not provide solution of equations (31) 
and (32) for the general case. However, we can give some 
recommendations and an example how these equations 
can be solved using method of characteristic function. 
We hope that with some modification this method can 
be used for other special cases as well. 



moment. Metallic spheres which also can be replaced by 
point dipoles are out of scope of our consideration since 
the orientation of their dipole moments depends on di- 
rection of external field. 



a) 




b) 



2r„ 
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III. LATTICE OF UNIAXIAL DIPOLAR 
SCATTERERS 

If the scatterers which form electromagnetic crystal 
are small as compared to the wavelength then sometimes 
they can be effectively replaced by point dipoles. It is 
assumed that the dipole moment of such a dipole is de- 
termined by local field acting to the scatterer and the 
field produced by the scatterer is equal to the field cre- 
ated by the dipole. The polarizability which relates the 
induced dipole moment and the local field acting to the 
scatterer is the only parameter which depends on the 
shape of scatterer in such a local field approach. 

Generally, the field produced by any scatterer can be 
presented using expansion by multipoles. The electric 
and magnetic dipoles are first and second order mul- 
tipoles. For some scatterers, the electric or magnetic 
dipole moments dominate over high-order multipoles. It 
means that some scatterers behave as electrical dipoles, 
some other ones as magnetic. Below we will consider 
only such scatterers. The scatterers which has both elec- 
tric and magnetic dipole moments of the same order or 
whose quadrupole or other high-order multipoles can not 
be neglected are out of scope of our consideration. More- 
over, below we will consider only scatterers which can be 
replaced by dipoles with fixed orientation. 

The typical example of the scatterer which behaves 
as electric dipole with fixed orientation at microwave fre- 
quencies is a short metallic cylinder or piece of wire which 
can be loaded by some inductance in order to increase its 
polarizability [15] (sec Fig. 4.b). At optical frequencies 
it can be prolate metallic cylinder which has strong plas- 
monic resonance. The typical magnetic scatterer at mi- 
crowave frequencies is split-ring resonator [16] (see Fig. 
4. a) if the bianisotropic properties of this scatterer are 
neglected or canceled using method suggested in [17]. At 
optical frequencies the split metallic rings [18, 19] behave 
as magnetic scatterers with fixed orientation of dipole 



FIG. 4: Geometries of scatterers which can be modeled as 
dipoles with fixed orientation: a) split-ring-resonator, b) in- 
ductively loaded wire. 



A. History of the problem 

Further in this section we present an analytical solu- 
tion for problem of plane wave diffraction on semi- infinite 
electromagnetic crystals formed by point scatterers with 
known polarizabilites, but before that we have to describe 
history of this problem. The first attempt to obtain such 
an analytical solution has been made by G.D. Mahan and 
G. Obermair in a seminal work [10]. Analytical expres- 
sions for reflection coefficients and amplitudes of excited 
modes for a semi-infinite crystal were obtained in terms 
of wave vectors of the infinite crystal eigenmodes. How- 
ever this theory is not free from the drawbacks. Mahan 
and Obermair treated the interaction between a refer- 
ence crystal plane of a semi-infinite crystal and its N 
nearest neighbors exactly, neglecting the other crystal 
planes. That is why this approximation is called "near- 
neighbor approximation" . Such an approach allows to 
introduce fictitious zero polarization at the imaginary 
crystal planes in free space over the semi-infinite crystal. 
This manipulation gave a set of equations which were 
treated in [10] as additional boundary conditions. It will 
be shown below that if the interaction between planes 
is taken into account exactly but not restricted to finite 
number of neighboring planes then the fictitious polar- 
ization of imaginary planes turns out to be non-zero. In 
works of C.A. Mead [20, 21] it was already shown that 
the "near-neighbor approximation" appears to be not a 
strict one. The mead states that the serious disagreement 
appears in the cases then the interaction between crys- 
tal planes falls off not sufficiently fast with distance. In 
other words, the results of Mahan and Obermair are valid 
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only then the high-order spatial Floquet harmonics pro- 
duced by the planes rapidly decay with distance. Mahan 
and Obermair considered only the normal incidence of 
the plane wave. Within such a restriction their approach 
is valid in the case when the periods of the structure 
are small as compared with the wavelength in the host 
medium. The strong disagreement with the exact solu- 
tion appears when one of high-order Floquet harmonic 
happens to be propagating one. This fact is illustrated 
below by a numerical comparison. 

The work [10] caused numerous extensions [22-24]. 
The Mahan and Obermair approach was generalized for 
the cases of oblique incidence [22], both possible polar- 
izations of incident wave [24] , various lattice structures of 
the crystal [22], tensorial polarizability of scatterers [24] 
and even diffraction of the finite-size slabs of the crystals 
was considered [23]. Note, that all the listed works use 
the same "near-neighbor approximation" and their ap- 
plicability is restricted as described above. In order to 
avoid this trouble one needs to use another model for in- 
teraction between crystal planes. The simplest one is the 
so-called "exp model" suggested by Mead [20, 25] which 
assumes that interaction can be described by a single de- 
caying exponent. In terms of spatial Floquet harmonics 
this approach is equivalent to neglecting all high-order 
Floquet harmonics except the one with the slowest decay. 
The "exp model" as well as the "near-neighbor approxi- 
mation" allow to solve the problem of excitation analyt- 
ically for both normal [20] and oblique [25] incidences. 
The "exp model" of Mead gives a set of two equations 
which correspond to the generalized Ewald-Oscen extinc- 
tion principle formulated in the present paper. The first 
equation of Mead is the same as one of equations given 
by "near-neighbor approximation" . It describes the fact 
that the incident electromagnetic wave (fundamental Flo- 
quet harmonic) inside the semi-infinite crystal is canceled 
by induced polarization of the crystal. This fact was 
pointed out in the papers [10, 20, 24]. The second equa- 
tion clearly expresses the fact that induced polarization 
cancels also the second Floquet harmonic (taken into ac- 
count in the "exp model") of incident wave which has 
zero amplitude, but unfortunately it was not noted by 
the authors. The system of these two equations is solved 
in [20] and the amplitudes of excited eigenmodes and an 
expression for reflection coefficient are obtained. 

It is possible to modify the "exp model" in order to 
obtain an exact solution. For that purpose one simply 
should take into account all Floquet harmonics in the in- 
teraction between the crystal planes. This has been done 
by authors of the present paper and the results are pre- 
sented below. As it was shown above, it turns out that 
every incident Floquet harmonic (even if it has zero am- 
plitude) is canceled by the induced polarization following 
to the generalized Ewald-Oseen extinction principle. It 
provides an infinite system of equations relating ampli- 
tudes of excited eigenmodes. This system can be trun- 
cated and then the number of equations in the system 
turns out to be equal to the number of Floquet harmon- 



ics taken into account. Such the finite system can be eas- 
ily solved analytically for the case when only two Floquet 
harmonics are taken into account (this is the "exp model" 
of Mead [20]), but in the case when one would like to take 
into account more Floquet harmonics this approach re- 
quires numerical calculations. We avoid the truncation of 
the system of equations and offer a closed-form rigorous 
analytical solution which is simple and explicit. 

Note, that a "formally closed solution" for the prob- 
lem under consideration was proposed by Mead in [21]. 
In this solution there is a contour integral of a certain 
function given in the form of infinite series. However, 
the calculation with help of such "formally closed solu- 
tion" requires serious numerical efforts. The main idea of 
work [21] is based on introduction of characteristic ana- 
lytical function which allows to determine all parameters 
entering the expression for the reflection coefficient. It is 
shown that knowledge of its roots allows to recover this 
function and obtain analytical expressions for all ampli- 
tudes of excited eigenmodes and for the reflection co- 
efficient, consequently. Unfortunately, these roots were 
not found in [21]. That is why the contour integration 
was used in [21] in order to bypass the problem of these 
roots finding. In fact, as it is shown below, the roots of 
this characteristic analytical function are determined by 
the wave vectors of Floquet harmonics and can be eas- 
ily expressed analytically. This fact is a consequence of 
generalized Ewald-Oseen extinction principle which is an 
important point of our theory. 

One could directly apply general results (26) and 
(28) to the problem under consideration. It is enough 
to replace polarization density of eigenmodes by three- 
dimensional delta function corresponding to the point 
of location of the dipolar scatterer in the unit cell and 
one could obtain the system of linear equations (26) and 
(28) which correspond to the boundary conditions in the 
present case as it is shown below. However, we prefer to 
re-deduce all the expressions by making the same steps 
as in the previous section but for the case of point scat- 
terers. We suppose that it is very useful step in order 
to demonstrate physical background of the computations 
undertaken in the previous section at a specific example. 



B. Dispersion equation 

Let us consider an infinite crystal formed by point di- 
electric dipoles with some known polarizability a along 
fixed direction given by unit vector d, a = add. The case 
of magnetic dipoles can be easily obtained from the the- 
ory for dielectric ones using duality principle. The scat- 
terers are arranged in the nodes of the three-dimensional 
lattice with an orthorhombic elementary cell a x b x c 
located in free space, see Figure 5. 

The distribution of dipole moments corresponding to 
an eigenmodc with wave vector q = (q x ,%,qz) T is de- 
scribed as p„ iSj z = p e ~i( q * an+q y bs+q * cl \ where n,s,l are 
integer indices of scatterers along the x— ,y— ,z— axes, 
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FIG. 5: Geometry of an infinite electromagnetic crystal 
formed by uniaxial dipolar scatterers 



respectively, and p is a dipole moment of the scatterer 
located at the center of coordinate system. Follow- 
ing to the local field approach p can be expressed as 
p = a(Ei oc d)d, where Ei oc is a local electric field acting 
to the scatterer. The local field is produced by all other 
scatterers which form the infinite crystal and it can be 
given by the formula 



E 



loc. — G(R IliSj ;)p ra . 



8,1 3 



(33) 



n,s,l 



where G(r) is the three-dimensional dyadic Green func- 
tion of the free space (9) and summation is taken over 
all triples of indices except the zero one. Accordingly 
the following dispersion equation for the crystal under 
consideration is obtained (compare with (6)): 



n,s,l 



G(R„, s ,z)e 



d. (34) 



In order to evaluate sums of series in (34) we use a 
plane- wise approach [26]. According to this approach 
the dispersion equation takes the following form: 



(35) 



The coefficients j3 n describe the interaction between 
planes and include the information on transverse wave 
vector components q y , q z as well as on a geometry of a 
single plane (lattice periods b, c). For n^O coefficients 
Pn can be expressed using expansion by Floquet harmon- 
ics. For n — the calculation of coefficient (3q (describing 
interaction inside of a plane and expressed in the form 
of two-dimensional series without the zero term) requires 
additional efforts (see [11, 26, 27] for details). 

The electric field produced by a single plane (namely 
two-dimensional grid b x c of point dipoles with the distri- 
bution p s j = pe~i( q y bs+qzCl, >) located in the plane x = 



is equal to 



E(r) = 



26c£ 



E < 

8,1 



slgn(x) 
I 



signO) 

8,1 



x pj_ 



where kf , = (±fc?„ 



(36) 

2irl 



K s — 1y+ b ' K l ~ 1 Z+ c ' 



a l = yfk 2 — (ks) 2 — (ki) 2 andfc is wave number of free 
space. One should choose the square root in the expres- 
sion for k s j so that Im(y^) < 0. The sign ± corresponds 
to half spaces x > and x < respectively. 

The formula (36) defines an expansion of the field pro- 
duced by a single grid of dipoles in terms of plane waves 
and it can be obtained using double Poisson summation 
formula to series of fields produced by single scatterers 
in free space. These plane waves have wave vectors . 
They are also called Floquet harmonics and represent a 
spatial spectrum of the field (compare with (14)). Flo- 
quet harmonics are widely used in analysis of phased ar- 
ray antennas [28]. 

Using (36) we get the following expression for (3 n (n =/= 
0): 



-sign(n) p —jk* it a\n\ 



(37) 



where 7^ = [k 2 - (k± ; • d) 2 )/ (2jbce k x s l ). After substitu- 
tion of (37) into (35), changing the order of summation 
and using the formula for sum of geometrical progression 
we obtain the dispersion equation in the following form: 



7 S 



7. j 



O i( fe ,.i+9x)o 



1 



i 



(38) 

This is a transcendental equation expressed by rapidly 
convergent series. Dispersion properties of the crystal 
under consideration can be studied with help of numer- 
ical solution of the latter equation. The case of dipoles 
oriented along one of the crystal axes of the orthorhom- 
bic crystal has been considered in [11], the dispersion 
equation was solved and typical dispersion curves and 
iso-frequency contours for resonant scatterers were pre- 
sented. 



C. Semi-infinite crystal 

Now let us consider a semi-infinite electromagnetic 
crystal, the half-space x > a filled by the crystal formed 
by point dipoles, see Figure 6. The structure is excited 
by a plane electromagnetic wave with the wave vector 
k = (k x , k y , k z ) T and the intensity of electric field Ei nc - 
Let us denote the component of the incident electric field 
along the direction of dipoles as i?i nc = (E inc • d). The 
axis x is assumed to be normal to the interface. The 
tangential (with respect to the interface) distribution of 
dipole moments in excited semi-infinite crystal is deter- 
mined by the tangential component of the incident wave 
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FIG. 6: Geometry of a semi-infinite electromagnetic crystal 
formed by uniaxial dipolar scatterers 



vector. It means that p n ,s,l — p n e^^ kybs+kzCl \ where 
the polarizations of zero-numbered scatterers from planes 
with the index n (parallel to the interface) are denoted 
as p n = p n ,o,o- The plane-to-plane distribution {p n } is 
unknown and it has to be found. Using the local field 
approach one can write the infinite linear system of equa- 
tions for this distribution: 



-jk x am 



ipn 



V m > 1. (39) 



D. Expansion by eigenmodes 

In order to solve (39) accurately one has to use an 
expansion of the polarization by eigenmodes [10]: 



(42) 



where A4 are amplitudes of eigenmodes and qi^ are the 
x-components of their wave vectors. Every eigenmode is 
assumed to be a solution of the dispersion equation (38) 

with the wave vector = (gi , k y , k z ) T . In the formula 
(42) the summation is taken by eigenmodes which cither 

transfer energy into half-space x > a (-Jjj- > 0) or decay 

along a:— axis (Im(gi^) < 0). 

Let us assume that the dispersion equation (38) is 
solved (for example numerically) and the necessary set of 

eigenmodes {q& } is found. Then the substitution of (42) 
into (39) will replace the set of unknown polarizations of 
planes by a set of unknown amplitudes of eigenmodes: 



-jq x am 



-^inc^ J^ am _j_ ^ ^ p n _ m S ^ Aie 



(43) 

Applying the auxiliary relation evidently following from 
(38): 



- E ^- 



the equation (43) can be transformed as follows: 



(44) 



The distribution {p n } of polarization in the excited semi- 
infinite crystal can be determined solving the system of 
equations (39). The known distribution of polarization 
allows to determine the scattered field in the half-space 
x < a with help of the expansion by Floquet harmonics 
(36): 



E 



E Es 



te 



(40) 



s,l 



where the amplitudes of Floquet harmonics are following: 



E, 



J 



+00 



2abeok 



— Kj x IK.i xd]] > p n e 



-jkg l an 



s.l 



n=l 



(41) 



If the crystal supports propagating modes, it is quite 
difficult to find a solution of (39) numerically. Simple 
methods such as system truncating (considering a slab 
with finite thickness instead of half-space like in [29]) 
results in nonconvergent oscillating solutions which have 
nothing to do with actual solution of (39). 



J2 A A E Pn- m e- iq * )an )=E iac e-^ am . (45) 



\n— — 00 



It should be noted, that using definition of Mahan and 
Obermair for the polarization of fictitious planes (p n — 



\M < 0) one can rewrite (45) as 



n=—c 



Pn—mPn -^inc^ 



-jk x am 



(46) 



It is evident that the assumption of Mahan and Ober- 
mair, requiring all polarizations of fictitious planes to be 
zeros, contradicts with (46). This fact proves that the 
"near-neighbor approximation" made in [10] is not accu- 
rate. 

The system of equations (45) can be truncated and 
solved numerically quite easily in contrast to (39). As a 
result, the amplitudes of eigenmodes {Ai} can be found 
and the polarization distribution can be restored using 
formula (42). The amplitudes of scattered Floquet har- 
monics (41) can be also expressed in terms of excited 
eigenmodes amplitudes by means of substitution of (42) 
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into (41), changing the order of summation and evaluat- 
ing sums of geometrical progressions. The final expres- 
sion for the amplitudes of scattered Floquet harmonics is 
following (compare with (28)) : 



E, 



Ki x [k-, x d]] 
2jabe k x l 



E* 



(47) 



One can stop at this stage and claim that the prob- 
lem of the semi-infinite electromagnetic crystal excita- 
tion is solved. However in this case the solution would 
require long numerical calculations, such as solving the 
system (45) and substituting the obtained solution into 
(47). The possibility to make all described operations 
analytically in the closed-form is shown below. 



E. Analytical solution 

Substituting the expansion (37) into (45) we obtain: 



eM e 



— E\ p -jk x am 
- tJ inc c ' 



|n— m| 



(48) 



Changing the order of summation in (48), taking into 
account that n — m < and using formula for the sum 
of geometrical progression we obtain: 



(49) 

This is a system of linear equations where unknowns are 
given by expressions in brackets. It has a unique solution 
because the determinant of the system has finite nonzero 
value. Note that k x 




k x only if (s,l) = (0,0). Thus, 



the solution of (49) has the following form: 



E 



A, 



1 



iW7o + o, if («, i) = (0, 0) 
0, if (s, 0^(0,0) 

(50) 

This equation could be directly obtained from general 
expression (26) by substitution of delta function instead 
of polarization density of eigenmodes, but as we already 
mentioned above we intentionally re-deduced it since we 
suppose that it can help to understand background for 
deduction of (26) for general case. The values at the right 
side of (50) are the normalized amplitudes of incident 
Floquet harmonics, and the series at the left side are the 
normalized amplitudes of Floquet harmonics produced 
by the whole semi-infinite crystal polarization which can- 
cel the incident harmonics. Thus, the equation (50) rep- 
resents the generalization of the Ewald-Oseen extinction 
principle already formulated above for general case of 



semi-infinite crystals: The polarization in a semi-infinite 
electromagnetic crystal excited by a plane wave is dis- 
tributed in such a way that it cancels the incident wave 
together with all high-order spatial harmonics associated 
with periodicity of the boundary. 

Note, that the formula for the amplitudes of scat- 
tered Floquet harmonics (47) contains series that have 
the same form as (50), but another sign in front of k* t . 

The amplitudes of the excited modes At can be found 
numerically from the infinite set of equations (50) and 
substitution of Ai into (47) will give us amplitudes of 
scattered Floquet harmonics. However, it is possible to 
obtain a closed- form analytical solution of (50). 

In order to solve the set of equations (50) one should 
consider a characteristic function f(x) (see also [21]) of 
the form : 



f(u) = uJ2 A i 



1 



u — eii* >a 



(51) 



Comparing (47) and (50) with (51) one can see that the 
function f(u) has the following properties: 

• It has poles at u = e J9 * >0 

• It has roots at u = e jh '-' a , (s, I) ^ (0, 0) and u = 

• It has a known value Ei nc /j^ at u = ei k * a 

• Its values at u = e~ jks - ia are equal to the normal- 
ized amplitudes of scattered Floquet harmonics 

• Its residues at u = e jq ^ >a are equal to the normal- 
ized amplitudes of excited eigenmodes. 

It is possible to restore the function f(u) using the 
known values of its poles, roots and a value at one point: 



/(«) = 



E inc u 



7oV J 



n 



JK,i a 



(s,0^(o,o) 



k x a, 



J k li a 



n 



— e jqx a 

(52) 

The knowledge of the characteristic function f(u) pro- 
vide us with complete solution of our diffraction problem. 
The amplitudes of excited eigenmodes with indices n are 
equal to residues of f(u) at u = e J9 ^" >a : 

-Einc(l 



An = Res/(w) 



■ (») 



7oto 



I [ _ Z _ Z TT 

II e jk x a _ e lK,i a 11 



(s,0#(0,0) 



(n) (i) ' 

< l ' a — e 3q * a 



(53) 



and the amplitudes of scattered Floquet harmonics with 
indices (r,t) can be expressed through values of f(u) at 



a = e 



-jk* t a. 



Er,f — 



E inc e-^ a [k~ t x [k~ x d]] 
2jabe k* t ^ fi ei k * a 
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n 



-oK,t a _ P iK,i a 



(«.0^(0,0) 



n 



-lK,t a 



(54) 



The products in the formulae (53) and (54) have very 
rapid convergence. It is enough to take a few terms 
in order to reach excellent accuracy. The main re- 
quirement for truncation of these infinite products is to 
take into account all terms corresponding to propagating 
Re(/Cg ;) = and slowly decaying Im(/cJ;) <C 2ir/a Flo- 

quet harmonics as well as propagating Re(gx ) = and 
slowly decaying Im(qx ) <^2ir/a eigenmodes. 



F. Comparison with other theories 

Let us consider the case from work [10] when d = 
yo, and a = b = c. In this case the formula (54) for 
the fundamental Floquct harmonic (r = t = 0) can be 
rewritten in terms of the reflection coefficient: 



R = - e ~ 2jk " a 



n 



-jk x a 



JK,l a 



(8,0^(o,o) 



e ] k x a 



JK i a 



n 



3 jk x a _ e jq^a 
-jk^a _ e jqi''a 



Comparing that result with the final result of the work 
[10] (the next formula after (C7) on page 841) one can 
see that the first product in our formula (55) 



IT 



n 



-jk x a _ jk* ,a 



(s,0^(0,0) 



(-:■> 



k x a 



JK i a 



(56) 



is absent in the result of Mahan and Obermair. This 
difference is a consequence of the fact that in our study 
we considered interaction between crystal planes accu- 
rately taking into account all Floquet harmonics for 
any distance between planes in contrast to the "near- 
neighbor approximation" used in the approach of Mahan 
and Obermair. 

The dependence of the product II vs. normalized fre- 
quency is plotted in Figure 7 for the case of normal in- 
cidence k y — k z = and k x — k. One can see that 
the value of the product is nearly equal to the unity for 
ka < 1.67T, but for ka > 2tt the value of the product 
significantly differs from the unity. Thus we conclude 
that the theory of Mahan and Obermair is valid in the 
low frequency range when periods of the lattice are small 
compared to the wavelength. Our theory does not have 
such a restriction (within the frame of the dipole model 
of electromagnetic crystal). 

The comparison with results of [21] shows that (55) 
is equivalent to formula (46) from [21] with II = exp(r) 
where T is given by the contour integral (47) from [21]. 
The calculation of II using (56) requires taking into ac- 
count only a few terms in the infinite products, because 
they are very rapidly convergent. This is a significant ad- 
vantage of our approach as compared to work [21] which 
requires complicate numerical calculation of the contour 
integral. 
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In the long- wavelength limit, the series in equation (34) 
for the cubic lattice can be replaced by the integral taken 
over the whole space except unit cell V and we obtain: 



a 



( 



G(R)e- j(qR) rfR d 



■ d. 



(57) 



The integral in the right-hand side of equation (57) can 
be evaluated by means of the same technique that was 
used while deducing Ewald-Oseen extinction principle in 
[9]. The result is following: 



(q-d) 2 -|q| : 
K 2 - Iql 2 



V 
so ' 



(58) 



The obtained dispersion equation (58) can be trans- 
formed in the common form: 



where 



£ ~(^-^)= £o (|q| 2 -^), 



a/(e V) \ 



so 1 



l-a/(3e V) J ' 



(59) 



(60) 



and qd = (q ■ d) is the component of the wave vector q 
along the anisotropy axis. 

The formula (59) is classical form of the dispersion 
equation for uniaxial dielectrics [9] with permittivity e 
along the anisotropy axis and e in the transverse plane. 
The expression (60) is the Clausius-Mossotti formula for 
the effective permittivity of cubic lattices of scatterers. 

In the long-wavelength limit the formula (54) for am- 
plitude of reflected wave simplifies as follows: 



Ef 



(E inc -d)[k- x [k-xd]] 

[k 2 - (k ■ d) 2 ] k x +q x 



(61) 
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where k = (—k x , k y , k z ) T is wave vector of reflected 
wave. The formula (61) represents a compact form of 
an expression for electric field amplitude of a wave re- 
flected from an interface between an isotropic dielectric 
and an uniaxial dielectric (see, e.g. [30]). Note, that 
in our case the situation is simplified as compared to 
the general case, because the incident wave comes from 
isotropic dielectric with permittivity e which is equal to 
the permittivity of uniaxial dielectric in transverse plane. 
It means, that an incident wave with normal polarization 
with respect to the anisotropy axis transforms at the in- 
terface in a refracted ordinary wave without reflection. 

Let us consider the reflection problem at the special 
case when d = yo, k y = and E; nc || yo- The nonzero 
components of the wave vector for the incident wave can 
be expressed in terms of incident angle 6i as k z = k sin 9i 
and k x — k cos From (59) we obtain that in this case 
the transmitted wave has x-component of wave vector 
equal to q x = ^ek 2 /e a — k 2 z = yje/e k cos 9 t , where 6 t 
in angle or refraction. With the result (61) we get the 
reflection coefficient in the form: 



R = 



k x — 1x ni cos 6i — ri2 cos 6t 



k x + q x Tii cos 0i + n 2 cos 9 t 



(62) 



where ni = ^/sojlo and ni = i/sjlo are indices of refrac- 
tion of the materials. The formula (62) coincide with the 
classical Fresnel equation [9] . This fact can be treated as 
an additional verification of presented theory. 



IV. LATTICE OF SPLIT- RING-RESONATORS 

In this section we apply theory presented in previous 
sections for study of reflection from semi-infinite cubic 
lattice of split-ring- resonators. 

The general dispersion equation (6) of the integral form 
in the case of point electric scattercrs transforms into 
transcendental equation (38). In the case when d = y 
the dispersion equation (38) can be rewritten in the fol- 
lowing closed form convenient for numerical calculations 
(see [11] for details): 



e a 1 (uj) = C(k,q,a,b,c), 



(63) 



where C(k, q, a, b, c) is dynamic interaction constant of 
the form: 



C(k, q, a, b, c) = - ^ ^ ^\ K ° ( PsC ^ cos (^ cn 



1=1 Re(p s )^0 



nb 



+ oo +OC 

+ E E 

s— — oo l — — oo 



-3 k s,l a 



cos q x a 



2jbck^ l cos ; a — cos q x a 



(64) 



Ro( Ps )= 



2/1 +°° 
Q 2bc\jkS fi ^_ 



c 

'Vi 



+ 



r s C J 
8tt 3 / 3 



47r6 3 



(2jkb + 3)s + 2- kbs 



■^ s 3 (s + 1)(s + 2 ) 



cos(q y bs) 



(jkb + 1) (4 logt+ + t 2 _ log*" + 2e 3kb cos(q y b)) 



-2jkb (t+ \ogt+ + t- log* - ) + (7 jkb + 3) 



and the following notations are used: 



Ps = \/(k y sf - k 2 , r s = 2q 2 z -pl 



-j(k+qz)c 



-j(k-q z )c 



oi{k-q z )c 



The calculations using (64) can be restricted to the 
real part of dynamic interaction constant C(fc,q, a, b, c) 
only, because its imaginary part is given by much simpler 
expression (see [11] for details): 



Im {C(k, q, a, b, c)} = 



fc 3 

67T 



(65) 



The series in (64) have excellent convergence that ensure 
very rapid numerical calculations. 

The case of magnetic scatterers can be considered us- 
ing duality principle. The expression (55) can be used 
for calculation of reflection coefficient by magnetic field 
(originally, this equation represented reflection coefficient 
by electric field). The dispersion equation (63) have to 
be rewritten for the case of magnetic point scatterers in 
the following form: 



HoaJiuj) = C(k,q,a,b,c), 



(66) 



where a m (u>) is magnetic polarizability of the scatterers. 
The analytical expressions for the magnetic polarizabil- 
ity a(oj) of split-ring-resonators with geometry plotted in 
Fig. 4 were derived and validated in [31]. The final result 
reads as follows [11]: 



a(u>) 



Ac 2 



(67) 



where A is amplitude, ojq is resonant frequency and 
r = Aujk 3 1 (6irfio) is the radiation reaction factor. The 
expressions for amplitude A and resonant frequency ujq 
in terms of dimensions of split-ring-resonators are avail- 
able in [11, 31]. In the present paper we will use typical 
parameters A — 0.1/xo« 3 and ujq = l/(a^/eoA i o)- 
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The dispersion properties of cubic lattice of split-ring- 
resonators with such parameters have been extensively 
studied in [11]. Using the theory of the present paper we 
will study reflection properties of such a metamaterial. 
Let us consider the case of cubic lattice (a = b = c) , nor- 
mal incidence (ky — k z — 0) and let the magnetic field 
of incident wave is along directions of magnetic dipoles 
d = yo. Numerical solution of dispersion equation (66) 
with q y = k y = and q z = k z = allows to get a set 
of wave vectors of excited eigenmodes {qf }. These wave 
vectors are plotted at the top of Fig. 8 as functions of nor- 
malized frequency ka. The point ka = 1 corresponds to 
the resonant frequency luq os split-ring-resonators. One 
can see that the propagating modes (Im^) = 0) exist 
only for ka < 0.978 and ka > 1.044. It means that a par- 
tial resonant band gap is observed for ka £ [0.978, 1.044]. 
At the frequencies inside of the band gap all the eigen- 
modes decay with distance. Note, that such decaying 
modes exist at all frequencies, not only inside of the band 
gap. The Fig. 8 shows only the eigenmodes with slowest 
decay llm^)! < l.Sir/a. There is an infinite number of 
other decaying modes which decay with distance more 
rapidly. The contribution of such modes into the reflec- 
tion coefficient is negligible as it was shown above. The 
decaying modes can be separated into the three classes: 

• evanescent modes, the modes which have Ke(q x ) = 
0, they decay exponentially from one crystal plane 
to the other one; 

• staggered modes, the modes which have Re(q x ) = 
Tt/a, they exponentially decay from one crystal 
plane to the other one by absolute value, but the 
dipoles in the neighboring plane are excited in out 
of phase; 

• complex modes, the most general case of the decay- 
ing modes which have He(q x ) ^ 0, they experience 
both exponential decay and phase variation from 
one crystal plane to the other one. 

The evanescent modes are the most common type of 
decaying modes. They can be observed in dielectrics with 
negative permittivity, for example. The staggered modes 
are limiting case of the complex modes and can be widely 
observed in periodical structures in vicinity of the band 
gap edges, see for example [32, 33]. The complex modes 
of general kind are quite exotic for common materials. 

In the system under consideration we are able to ob- 
serve all three kind of mentioned decaying modes. The 
presence of staggered and complex modes are evidence of 
spatial dispersion in this material reported in [11]. The 
staggered modes exists for ka > 0.984, evanescent modes 
for ka > 1.015 and complex modes for ka S [0.984, 1.015] 
(see Fig. 8). One can see that for a fixed frequency 
from the range ka £ [0.978, 0.984] there are two stag- 
gered modes and in the range ka £ [1.015,1.044] there 
are two evanescent modes. Actually, in the range ka € 
[0.984, 1.015] there are also two complex modes which 
have the same imaginary parts but differs by sign of the 
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FIG. 8: (Color online) Dependencies of the normalized wave 
vectors q x a/n of excited eigenmodes (imaginary and real 
parts) and the reflection coefficient R calculated using (55) 
on normalized frequency ka for semi-infinite cubic lattice of 
split-ring-resonators with A — 0.1/ioa 3 and loq = l/(a v /£ °A t o)- 



real part. Thus, we conclude that for every frequency 
from the range ka G [0.95, 1.08] that we consider the in- 
cident wave will excite in the crystal a pair of modes 
with |Im(g.i;)| < 1.57r/a (propagating and staggered, two 
staggered, two complex, two evanescent or propagating 
and evanescent). Using the usual approach one has to 
introduce an additional boundary condition to solve this 
problem since usual condition of tangential component 
continuity is not enough in the case of excitation of two 
modes. 

Using the theory introduced in the previous section it 
is enough to substitute obtained wave vectors of eigen- 
modes into (55) the reflection coefficient is calculated. 
The reflection coefficient is plotted at the bottom of Fig. 
8. One can see that at the frequencies in vicinity of the 
bottom and top edges of the band gap the semi-infinite 
crystal operate nearly as the electric and magnetic walls, 
respectively, as it was predicted in [32] . At the frequency 
ka = 0.984 the reflection coefficient is equal to +1 (elec- 
tric wall), and at ka — 1.044 it is —0.8 + 0.6j (nearly 
magnetic wall). Note, that the frequency correspond- 
ing to the electric wall effect is not equal to the bot- 
tom edge of the band gap and there are no frequency 
exactly corresponding to the magnetic wall effect. The 
use of the usual formulae for reflection coefficient from 
magnetic and Clausius-Mossotti formulae which do not 
take into account effects of spatial dispersion one could 
get idea that magnetic and magnetic wall effects have to 
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happen at the edges of the band gap. Our study demon- 
strate that if the spatial dispersion is taken into account 
accurately then it is not so. 

Thus, we have demonstrated how the proposed the- 
ory can be used for modeling of reflection from semi- 
infinite crystals with spatial dispersion. Our theory can 
be treated as generalization of results of Mahan and 
Obermair [10] which have been widely applied for model- 
ing of various kinds of reflection problems. We hope that 
the present generalization can find much more applica- 
tions in modeling of reflection from spatially dispersive 
materials since in has no restriction on the period of the 
lattice to be smaller than wavelength and allows to con- 
sider electromagnetic crsytals of general kind. 

V. CONCLUSION 

In this paper a new approach for solving problems 
of plane- wave diffraction on semi-infinite electromagnetic 
crystals is proposed. The boundary conditions for the in- 
terface between isotropic dielectric and electromagnetic 
crystal of general kind are deduced in the form of infi- 
nite system of equations relating amplitudes of incident 
wave, excited eigenmodes and scattered spatial harmon- 
ics. This system of equations represent mathematical 
content of generalized Ewald-Oseen extinction principle 
which is formulated in this paper: the polarization of the 
semi-crystal excited by plane wave is distributed in such 
way that it cancels the incident wave together with all 
high-order spatial harmonics associated with periodicity 
of the boundary. In our opinion, the proof of general- 



ized Ewald-Oseen extinction principle presented in this 
paper is an important theoretical fact which helps to un- 
derstand interrelation between reflection and dispersion 
properties of electromagnetic crystals. If the eigenmodes 
of the infinite crystal are known then the system can be 
solved numerically which provides new numerical method 
for solving the diffraction problem under consideration. 
We believe in the quite good prospects for the applica- 
tion of the described method in further studies of dielec- 
tric and even metallic electromagnetic crystals at both 
microwave and optical ranges. 

For the special case when the crystal is formed by small 
scatterers which can be effectively replaced by dipolcs 
with fixed orientation the deduced system of equations 
is solved analytically using method of the characteris- 
tic function. The closed form expressions for the ampli- 
tudes of excited eigenmodes and scattered spatial har- 
monics are provided in terms of rapidly convergent prod- 
ucts. These expressions can be treated as generaliza- 
tion of classical result of Mahan and Obermair [10] for 
the case when period of the lattice can be large as com- 
pared to the wavelength. The proposed method is ap- 
plied for calculation of reflection coefficient from semi- 
infinite crystal formed by resonant magnetic scatterers 
(split-ring-rcsonators) at the frequencies corresponding 
to the strong spatial dispersion. 
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